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In this paper, the linear Gaussian relay problem is considered. Under the linear time-invariant (LTI) model 
the problem is formulated in the frequency domain based on the Toeplitz distribution theorem. Under the further 

00 

■ assumption of realizable input spectra, the LTI Gaussian relay problem is converted to a joint design problem of 

source and relay filters under two power constraints: one at the source and the other at the relay, and a practical 
solution to this problem is proposed based on the projected subgradient method. Numerical results show that the 



proposed method yields a noticeable gain over the instantaneous amplify-and-forward (AF) scheme in inter-symbol 
interference (ISI) channels. Alsc 
established in flat-fading channels 



interference (ISI) channels. Also, the optimality of the AF scheme within the class of one-tap relay filters is 
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I. Introduction 

Relay networks have drawn extensive interest from research communities because they play an im- 
portant role in enlarging the network coverage in wireless communications. Although the capacity of 
relay networks is not exactly known yet, many ingenious coding strategies including decode-and-forward 
(DF), compress-and-forward (CF), etc. beyond simple AF schemes have been developed [1,2]. Recently, 
Zahedi et al. proposed an advanced linear scheme for relay networks based on (strictly-)causal linear 
processing at the relay to compromise the complexity and performance between the complicated coding 
strategies and the simple AF* scheme [3,4]. While information theorists approached the problem from 
the perspective of capacity and capacity-achieving schemes [5-8], researchers in the signal -processing 
community also tackled this problem based on measures like the received signal-to-noise ratio (SNR) or 
minimum mean square error (MMSE). Most of their results are based on the setup in which linear pro- 
cessing is at the relay and destination but not at the source, e.g., [9, 10]. Although these works provide 
meaningful approaches to the relay problem, it is not optimal not to have processing at the source from 
the fundamental perspective of data-rate maximization. To this end the processing at the source such 
as the input covariance function design should be incorporated together with the processing at the relay. 
(Once the processing at the source and relay is fixed, the optimal destination processing is automatically 
given for several well-known criteria.) However, the joint design of source and relay processing is a hard 
problem even in the linear Gaussian case, as shown in [3,4]. In [3,4], the authors considered general 
time-varying linear processing at the relay in Gaussian channels. Although they obtained the capacity 
for frequency-division strictly-causal linear relaying, the general linear relay case was not explored fully 
[3,4]. In the general linear relay case, the problem is a sequence of non-convex optimization problems, 
and it is seemingly intractable. To circumvent such difficulty, in this paper we consider tractable and 
practical LTI filtering at the source and relay. We find that it is still a hard problem to obtain the capacity 
with a single-letter characterization even in this case because the search space still has countably infinite 
dimensions; optimal source and relay filters may have infinite impulse responses (IIRs). However, we 
provide a practical solution to design the source and relay filters jointly to maximize the transmission 
rate for general ISI Gaussian relay networks. 

Under the LTI framework, the linear Gaussian relay problem can be formulated in the frequency do- 
main using the Toeplitz distribution theorem [1 1, 12]. When the relay filter is given and there is no power 

*In this paper, the AF scheme means the instantaneous AF scheme, which can easily be implemented by simple analog 
processing. 
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constraint on the relay, the problem reduces to the classical ISI channel problem for which the optimal 
strategy is known as water-filling in the frequency domain [13, pp. 407 - 430]. However, the freedom 
to design the relay filter and the power constraint at the relay make the problem far more difficult than 
the classical ISI channel problem, especially when stability and causality constraints are imposed on the 
source and relay filters. Our approach to this problem is that we first convert the problem to a constrained 
optimization problem in a finite dimensional space by restricting the source and relay filters to the class of 
finite impulse response (FIR) filters as in most practical filtering applications, and then apply the projected 
subgradient method, initially proposed by Polyak [14] and fully developed by Yamada et al. [15, 16], to 
this problem. Numerical results show that our method performs well and yields a noticeable gain over 
the AF scheme in ISI relay channels. 

Notations and Organizations 

We will make use of standard notational conventions. Vectors and matrices are written in boldface with 
matrices in capitals. All vectors are column vectors. For a scalar a, a* denotes its complex conjugate. For 
a matrix A, A T , A^ and tr(A) indicate the transpose, Hermitian transpose and trace of A, respectively, 
and A(m, n) denotes the m-th row and n-th column element of A. diag(di, • • • , d n ) denotes a diagonal 
matrix with elements d\, ■ ■ ■ ,d n . I n stands for the identity matrix of size n (the subscript is omitted 
when unnecessary), and denotes a vector of all zero elements. For a vector a, ||a|| denotes its 2-norm. 
The notation x ~ S) means that x is Gaussian-distributed with mean vector /x and covariance 

matrix S. E{-} denotes the expectation. For two signal processes x[n] and y[n], x[n] * y[n] denotes the 
convolution of the two processes. IR, I, I + and N denote the sets of real numbers, integers, nonnegative 
integers and natural numbers, respectively. For two sets A and B, A\B denotes the set minus operation. 



This paper is organized as follows. The system model and background are described in Section II. In 
Section III, the rate formula in the frequency domain is derived under the LTI model, and the performance 
of LTI relaying in flat-fading channels is investigated in Section IV. In Section V, a joint source and relay 
filter design method is proposed based on the projected subgradient method, and its performance in ISI 
channels is examined in Section VI, followed by conclusions in Section VII. 



We consider the general discrete-time additive white Gaussian noise (AWGN) relay network composed 
of source, relay and destination nodes, as shown in Fig. 1, where the source and relay nodes have maxi- 
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II. System Model and Background 
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mum available average power P s and P r , respectively. We assume that all propagation channels (i.e., the 
source-to-relay (S-R), relay-to-destination (R-D) and source-to-destination (S-D) channels) are linear, 
time-invariant and causal, and their impulse responses are absolutely summable, i.e., l^sr[^]| < oo, 
< 00 an d S£ol^sd[(]| < co> where h sr [l], h rd [l] and h sd {l] are the S-R, R-D and 
S-D channel impulse responses, respectively. Due to the absolute summability, the z-transforms of 
the propagation channel impulse responses are well-defined and given by H sr (z) = X^£o hsr[l]z~ l , 
H r d{z) = YlhLohrdl^z' 1 and H sd {z) = Y^h=Q^sd[l]z~ l ■ Then, the received signals at the relay and 
destination at the n-th symbol time are given by 

y r [n] = h sr [n] * x s [n] + w r [n] , and (1) 
Vd[n] = h sd [n] * x s [n] + h rd [n] * x r [n] + w d [n], (2) 

respectively, where x s [n] is the transmitted signal process at the source; x r [n] and y r [n] are the trans- 
mitted and received signal processes at the relay, respectively; y d [n] is the received signal process at 
the destination; and the noise processes w r [n] at the relay and Wd [n] at the destination are independent 
zero-mean white Gaussian processes with variance a 2 . 



w r [n] 




Fig. 1 
System model 



We consider the linear and causal processing at the relay. The general causal linear processing at the 
relay is given by 

x r [n] = ^dniyrll], (3) 

Kn 

for arbitrary linear combination coefficients d n i, as considered in [3,4]. However, such linear processing 
requires time-varying filtering at the relay, and is not readily realizable. Thus, in this paper, we restrict 
ourselves to the case of LTI causal filtering at the relay, as shown in Fig. 2. In this case, the relay output 
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«VN Relay 
Fig. 2 

System model with linear time-invariant filters 



is given by 



x r[n] = y^MVr[ n ~ J]> 



(4) 



1=0 



where [ho , hi , /t2 , ■ ■ ■ ] is the time-invariant impulse response of the relay filter and its z-transform is given 
by H(z) = Ylh^o hi z ~ l - ( m trie case of strict causality, we have ho = 0.) The received signal (4) at the 
relay can be written in matrix form as (5), and the filtering matrix H n in (5) has a Toeplitz structure. 



x r [0] 
x r [l] 

x r [n — 1] 



ho 

hi h ••• 

hi h\ ho : 



ho hi hn 







w r [0] 


y r [i] 
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_ Vr[n - 1] _ 




w r [n — 1} 



(5) 



:H ? , 



We assume the stability (i.e., X^£o 1^*1 ^ 00 ) an( ^ realizability^ for the relay filter. Since all processing 
from the source and to the destination is linear and time-invariant, the received signal at the destination 
in the z-domain is given by 



Y d (z) = (H sd (z) + H rd {z)H{z)H sr {z))X s (z) + H rd {z)H{z)W r {z) + W d (z), 



(6) 



where W r (z) and W d (z) are the z-transforms of noise processes w r [n] and w d [n], respectively. 

^It means that the LTI response H(z) has a rational transfer function and it can readily be implemented by an autoregressive 
moving-average (ARM A) filter [17]. 
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A. Background 

In this subsection, we briefly summarize some relevant results including the eigen-structure of Toeplitz 
matrices and the spectral factorization for the development in later sections. For a zero-mean* stationary 
random process y[n], the covariance sequence and its z-spectrum are given by 

oo 

r y [k) = E{y[n]y*[n - k}} = r* y [-k] and S y {z) := £ r y [k\z~ k = E{Y(z)y[0}}, (7) 

k=— oo 

respectively, where Y(z) = ^ n y[n]z~ n . The covariance matrix of a finite collection y n := [y[0], y[l], • • • 
y[n — 1]] T is given by 



Y.I :=E{y n y^} 



r„[0] r y [-l] 
r y [l] r y [0] 



r y [n-l] r y [n-2] 



r y [-n+l] 



r y [-l] 



r y [0] 



(8) 



Theorem 1 (Asymptotic eigen-structure of Toeplitz covariance matrices [12], p. 135) Let r y [k] be an 
absolutely summable autocovariance sequence of a stationary process y[n], let S y {e^) be its power 
spectral density (PSD), i.e., S y (e JU ) = S y ( z)\ z=e ju , and let D n be the n X n matrix, 

D n = diag(5 y {e>°),S v (e^ 1 ) , S y (e~^ ),■■■, S y (e^~^ ) , S y (e~^-^ ) ) , 

where uj\ = — . Then, for the covariance matrix H y n = [r y [k — l]]^ l=1 , the components of X^ n ) := 
WS-[W _1 — D n converge to zero uniformly as n — > oo (i.e. sup 1<fc /<„ |X( n )(fc, l)\ — > 0), where W 
is the discrete Fourier transform (DFT) matrix. 

For even n, we have a similar result with a slight modification. Theorem 1 simply states that the eigen- 
values of the Toeplitz covariance matrix of a stationary process are the uniform samples of its spectrum. 
Using Theorem 1, the following can easily be shown. 

Theorem 2 (Toeplitz distribution theorem [1 1], p. 65) Let {A,- r ^} be the eigenvalues of the Toeplitz co- 
variance matrix of a stationary process y[n\. Then, 



lim -j2f(X ( r ) ) = ^ T f(S : 
i=l J ~ n 



for any continuous function /(•). 
* In the case of a known non-zero mean, the mean of the process can be subtracted and the result can still be applied. 



(9) 
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In addition to the asymptotic eigen-structure of Toeplitz covariance matrices, we need some background 
in the spectral theory for stationary random processes, especially canonical spectral factorization. 

Definition 1 (Canonical Spectral Factorization [18], p. 197) Let S y (z) be a rational z-spectrum of a 
finite power process and assume that S y (z) is strictly positive. Then, the canonical spectral factorization 
of S y (z) is given by 

S y {z) = L{z) le L\z), (10) 

where L(z) = * s a un iq ue stable, causal, monic and minimum-phase (SCAMP) filter (i.e., the 

zeros and poles of L(z) are strictly inside the unit circle and L(oo) = 1 (or equivalently Iq = 1)), and 
7 e > 0. Here, L$(z) := L*(z~*) denotes the para-Hermitian conjugate. 

III. The Rate Formula in Frequency-Domain for LTI Relays 

First, note that the overall channel model (6) with LTI relay filtering is still a linear additive stationary 
Gaussian noise channel. Thus, for a given relay filter, the overall channel with the LTI relay filter reduces 
back to the classical ISI channel with stationary Gaussian noised In this case, stationary Gaussian signal 
processes with well-defined spectra are sufficient to achieve the capacity [13, pp. 407 - 430]. Hence, 
we assume that the source (or input) process x s [n] is a stationary Gaussian process. By concatenating 
symbols at the source up to time n — 1, we have 

<:=[x a [0],x a [l],... ,x s [n-l]f ~JV(0,E*«), C 11 ) 

and vectors y r n , and y^ are constructed similarly for the relay and destination nodes. Then, the power 
constraints for the source and relay are respectively given by 

(l/n)tr(S*«) < P., and (12) 
(l/n)E{tr(H n y"(H n y") H )} = (l/n)tr(H n (H^S^H^ + <r 2 I)Hf ) < P r , (13) 

where H* r is the filtering matrix for the S-R channel constructed based on {/i sr [Z]} similar to H n in (5). 
Thus, the maximum rate with LTI relaying for block size n is given by maximizing the mutual information 
between x* and y^ over and H n under power constraints (12) and (13), and the capacity with LTI 
relaying is given by its limit 

C LT i= lim sup -J(x*;y£), (14) 

n->oo sr^xs n 

§ However, the major difference between the two problems is that in the relay problem we even have to design the overall 
channel by properly choosing the relay filter and the power constraints at the source and relay are interwined. 
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sr \H I 



= log |a 2 H; d H n H^ (H^ + <r 2 I + (H* n d + H^H n H-)E^(H^ + H^H n H-) J 

- log \a 2 W n d H n H* (H r n d ) H + a 2 I\. (15) 

Here, Hff and are the filtering matrices for the S-D and R-D channels, respectively. Note that (14) is 
still valid for general linear time- varying relay filtering with H n given by an arbitrary lower triangular ma- 
trix. As mentioned in [4], the computation of capacity and the design of capacity-achieving (or at least rea- 
sonable) and H n are difficult problems in the case of general linear causal relay filtering. In the time- 
varying case, if we increase n by one, at least 2n new variables (n, 1), (n, 2), • • • , 5]^ s (n, n), 
dni, d n 2, • • • , d nn } appear (see (3)), and thus the complexity of the problem increases with the order of 
n! to make the problem difficult [3,4]. In the LTI case with a stationary source process, however, we have 
only two new variables r Xa [n — 1] and /i n _i for the increase of the problem size from n — 1 to n because 
of the Toeplitz structure of the covariance matrix in (8) and the filtering matrix in (5). Following the best 
input covariance matrix and relay filter for the problem size n is equivalent to designing the best infinitely 
long autocovariance sequence {r^ [A;] , k = 0, 1, ■ ■ ■ } and infinitely long relay filter {hi, I = 0, 1, ■ ■ ■ } 
first and then increasing the problem size. Thus, in the LTI case, we have 

1 r 



Clti = sup lim 

{r Xs [k]},H(z) n ^ n 



I(K-y n )\K°(b x Ak]}),H n (H(z)) 



(16) 



where the respective dependence of S™ s and H n on {r Xa [k]} and H (z) is explicitly shown. Here, taking 
the limit of n simplifies the problem significantly due to Theorem 2 since the eigenvalues are strictly 
positive due to the additive noise term and since f(t) = log(i) is a continuous function of t for t > 0. By 
Theorems 1 and 2 we have 

C LT i = sup — / - log 2 1 + S Xa (e 3 ) du, 

(17) 

where the input spectrum S Xa {e^) = Y^k=-oo r ^ B [k]e~^ k , since the eigenvalues of a covariance matrix 
are the samples of its spectrum and the determinant of a covariance matrix is the product of its eigenvalues. 
Here, we define the overall channel-to-noise power ratio (CNR) density as 

[€ h <r 2 (\H rd (ei")H(ei")\ 2 + l) D(e^) : K ' 

where N(e^ UJ ) and D{e^) are the numerator and denominator of the CNR density, respectively. Note that 
the CNR density captures the overall channel response from source to destination. When the CNR density 
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is multiplied by the input signal PSD, the product becomes the overall SNR density at the destination. 
(This quantity will be used in later sections.) In addition to the rate formula (17) in the frequency domain, 
the power constraints can also be expressed in the frequency domain as n — > oo. As n — > oo, again by 
Theorems 1 and 2, the power constraints (12) and (13) are respectively given by 

^ J S Xa (e j ")dw < P s , and (19) 

±J* \H(en\ 2 (\Hsr(en\ 2 S Xs (e juJ ) + v 2 )dto < P r , (20) 

since the trace of a matrix is the sum of its eigenvalues. Thus, the LTI relay problem is summarized by 
(17), (19) and (20). Note that for a given relay filter H(z) the problem without the power constraint (20) 
reduces to the well-known ISI channel problem and the solution of S Xa {e^ w ) is given by water-filling in 
the frequency domain [19]. However, the freedom to design H(z) and the relay power constraint (20) 
make the problem far more difficult than the simple ISI channel problem. To construct a practical method 
to solve this problem, we further assume that the input spectrum S Xa (e^) is also realizable. That is, its 
canonical spectral factorization is given by 

S Xs (z) = af(z)f^z) = T(z)T^z), (T(z) = y/Zr(z)), (21) 

where the SCAMP filter T(z) has a rational transfer function and, thus, S Xs (z) is a rational spectrum. 
In this case, the source process x s [n] can be modelled as the output of the stable and causal ARMA 
filter T(z) driven by a white Gaussian process x s [n] with unit variance, as seen in Fig. 2. Thus, the 
rate maximization problem under LTI relaying with realizable input spectra now reduces to a. joint design 
problem of LTI source and relay filters. Obtaining the capacity in a closed form still seems to be a difficult 
problem even in the LTI relay case. However, we propose a very effective and practical solution to this 
joint filter design problem in Section V. Before we tackle this problem, we investigate the problem in the 
case that all S-D, S-R and R-D channels have flat frequency responses in the next section. 

IV. Examination of LTI Relaying in Flat-Fading Channels 

In the case of flat fading, we have the system model (6) in which each of S-R, R-D and S-D channels 
has only one tap, i.e., H s d(z) = 1, H sr (z) = a and H rd {z) = b, as considered in [3,4]. Then, the 
received signal model in the z-domain is given by 

Y d {z) = (1 + abH{z))X s (z) + bH(z)W r {z) + W d (z). (22) 
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A. The One -Tap Relay Filter Case 

First, consider the well-known AF relaying. In this case, we have 

x r [n] = dy r [n], 



where E{x^} = P s and < d < J a ip^ +a i ='■ d max to satisfy the power constraints, and the received 
signal model is given by 

yd[n] = (1 + abd)x s [n] + bdw r [n] + (23) 

Due to the simple data model (23), the achievable rate in this case is known and given by 

1 / (l + abd) 2 P s \ 

R AF = max -log 1+ V ; ± , 24 

o<d<d max 2 \ b 2 d 2 + l a 2 ) 

and the optimal value of d is explicitly given by 



Vs^b}- (25) 

Now consider the one-tap LTI relay filter with an arbitrary delay: 

H(z) = dz~ A (26) 

for some integer A > 0^, where the relay gain d can be optimized under the power constraints. Note in 
the system model (22) that the relay filter affects both the channel gain and noise spectrum. However, in 
the one-tap relay filter case, the problem is simplified because the overall noise spectrum is white. In this 
case, the overall channel gain is given by 1 + abH(z) = 1 + abdz~ A and the overall noise spectrum is 
given by 

b 2 H(z)H\z)a 2 + a 2 = 6V 2 (^" A )(^- A )« + a 2 , 
= (b 2 d 2 + l)a 2 , 

since (z~ A )^ = z A . Note that the overall noise process in this case is white and equivalent to that in the 
AF data model (23); both have the same variance (b 2 d 2 + l)a 2 . Thus, the spectrum of Yd(z) is given by 

S Vd (e> u ) = \1 + abH(en\ 2 S Xs (en + (b 2 d 2 + l)a 2 , (27) 
and the channel frequency response is explicitly given by a raised-cosine function: 

\l + abH(e jul )\ 2 = (l + abde- ju,A )(l + abde julA ), 

= 1 + 2abdcos(ojA) + a 2 b 2 d 2 > 0. (28) 

^Note that strict causality implies A > 0. 
August 9, 2011 DRAFT 
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Since \H{e^)\ 2 = d 2 for H{z) = dz~ A , from (19) and (20) the power constraints are given by 



2tt 

J2 



-!- / S Xa {e? u )du <P S , and (29) 



|- ^ (a 2 5 Xs (e^) + a 2 )du = d 2 (a 2 ±- j* S Xs (e?")du> + a 2 ^j < P r , (30) 

which are the same as those of the AF scheme with A = 0. 

The problem with the given relay filter H(z) = dz~ A reduces to the simple ISI channel problem, and 
the optimal input spectrum S* g (e^) is obtained by water-filling under the two simple power constraints 
(29) and (30). In the following theorem, we establish the optimality of the AF scheme within the class of 
all one-tap relay filters. 

Theorem 3: Among all one-tap linear relay filters, i.e., H(z) = dz~ A with A € I, the AF scheme with 
A = maximizes the achievable rate. 

Proof: For a given A / 0, let 

VU*),* ) = argmax - y log (l + 1 ; V + 1) J S„W")) to. (3D 



Then, 



< ht^( 1 + W^^^ )]du ' (32) 



= Raf- (35) 

Here, (32) is obtained because |1 + abd*e~^ A \ 2 < (1 + abd*) 2 . (33) is obtained because the feasible 
set (S Xs (e JUJ ),d) satisfying the power constraint for A ^ is the same as that for A = when H(z) = 
dz~ A . (See (29) and (30).) (34) is obtained by Jensen's inequality. Finally, (35) is obtained by the 
definition of Raf in (24). ■ 

Theorem 3 states that the AF scheme with A = performs best within the class of one-tap relay filters 
with arbitrary delays. This is because the AF scheme achieves coherent signal combining between the two 
signal paths S-D and S-R-D. Instead of using the optimal water-filling source filter, we can also consider 
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a simple channel-equalizing source filter. However, the performance in this case is bad, as shown in the 
following theorem. 

Theorem 4: The achievable rate by an equalizing source filter for the one-tap relay filter H(z) = dz~ A 
is given by 

1 / \-(abd) 2 P s \ 

R\-tap,EQ = SUP - log \l + -p^j- • ~ 2 ) 06) 

regardless of the value of A > 0. Further, the supremum is given by Ri-tap,EQ = \ log (l + ^i) 
achieved when d = 0. 

Proof: We have X s (z) = yfaT{z)X s (z), where X s {z) is the z-transform of the white Gaussian process 
x s [n] with unit variance, and the equalizing source filter is given by 

f « = IT^im = TT^F* = 1 - oMz " + <oM)2z " 2A - ( " M)3z " 3i + • • • • 

When < d < X, the overall channel response 1 + abH(z) is SCAMP and, thus, the channel-equalizing 
source filter T(z) is also SCAMP. By the power constraint at the source, we have 

i r a r i 

— / S x .(e? u )du> = — — — T^2 ^ = P s ( 37 ) 

Z7r J_ n 2tt J_ n 1 + 2abdcos(ujA) + (abd) z 

because S x ,(z) = af(z)f^z). Since 1+2aMcos( L A )+(afcd)^ a; = I^gbp for every integer A > 0, 
a = (1 — (abd) 2 )P s regardless of the value of A > 0. With the channel-equalizing source filter T(z) = 
y/aT(z), the data model is given by yd[n] = yfax s [n\ + w e ff[n), where x s [n) ~ A/"(0, 1) and zu e // [n] ~ 
A/ r (0, (6 2 d 2 + l)cr 2 ), and the con - esponding achievable rate is given by (36). Now consider CNR(d) = 
^ffi+r in (36). Its derivative with respect to (w.r.t.) d is given by CNR'(d) = \-2a 2 b 2 d(b 2 d 2 + 1) - 
(1 - a 2 b 2 d 2 )2b 2 d]/[(b 2 d 2 + l) 2 cr 2 ] < for all d > 0. Thus, the rate is maximized when d = 0. ■ 
Theorem 4 states that it is optimal to turn off the (one-tap) relay filter when the channel-equalizing filter 
is to be used at the source. Thus, using the channel-equalizing source filter is not a proper choice for relay 
networks. 

Fig. 3 shows the achievable rates of several relaying schemes. For Fig. 3 (a) and (b), which show the 
same curves with two different x-axis ranges, we set a = 1 and b = 2, as in [4]. It is seen that simple linear 

To 1 

strictly causal schemes (one based on the filtering matrix in [4] and the other based on one-tap 

d 

filtering H(z) = dz~ l ) can outperform the CF scheme in the low SNR region, as already known from 
[4]. In this case of a = 1 and 6 = 2, the AF scheme achieves the cut-set upperbound f or P s = P r > 1 /3 
[6, Proposition 9]. It is interesting to observe that the simple linear scheme in [4] with n = 2 performs 
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Fig. 3 

Achievable rates of several schemes in flat-fading channels (P s = P r , a 2 = 1) - (a) and (b): 

a = 1, b = 2, AND (C) AND (D): a = 2, b = 2 



better than H(z) = dz^ 1 filtering in some low SNR values, although the latter outperforms the former 
eventually at high SNR. Fig. 3 (c) and (d), again showing the same curves in two different x-axis ranges, 
show the achievable rates when a = b = 2. In this case, it is seen that there is a gap between the cut-set 
bound and the AF scheme. In all the cases, it is seen that the two strictly causal linear schemes (one based 
on two-symbol concatenation in [4] and the other based on one-tap LTI filtering H(z) = dz^ 1 ) do not 
outperform the AF scheme, as expected by Theorem 3. 
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B. The Multiple-Tap Relay Filter Case: Insights from Ideal Low-Pass Filtering Relays 

In Section IV-A, it is shown that one-tap relay niters do not outperform the AF scheme in flat-fading 
channels. This is because any one-tap relay filter with a causal or non-causal non-zero delay cannot 
change the noise spectrum, but destroys the coherent signal combining that is available in the AF scheme. 
However, this is not the case when the relay filter has multiple taps. In this case, the overall noise spectrum 
as well as the channel gain spectrum in (22) can be shaped by the relay filter, and the LTI relaying scheme 
with multiple taps can outperform the AF scheme in flat-fading channels. However, the performance 
analysis in this case is far more difficult than that in the one-tap relay case, especially when the causality 
constraint is imposed on the relay filter. To circumvent this difficulty, in this subsection we relax the 
causality constraint on H(z) and consider the tractable ideal low-pass" relay filtering to gain insights into 
the interaction between the source and relay filters and to assess the rate gain that can be obtained by 
multiple-tap relay filtering. The frequency response of the ideal low-pass relay filter is given by 



H(e ju} ) 



M < 0J r 



0, UJ C < \U)\ < IT, 



(38) 



where 6 is the passband gain and uj c is the cutoff frequency. For a given to c , the optimization problem 
(17,19,20) with the ideal low-pass relay filter is expressed as 



max 

<5,S Xa (eJ-) 



subject to 



(1 + ab8) 2 a . , \ , 1 f 1 



1 



(39) 



1 

7T Jo 



S x .(e? u )dbJ - P. < 

T -/o 

5 2 {a 2 S Xs {e juJ ) + a 2 )dw - P r < 

S Xa (e jul ) > 0, VwG [0,tt], 



(40) 

(41) 
(42) 



where the even symmetry of spectra is used. Note that the problem is not jointly convex w.r.t. 5 and 
Sxsi^^) for a given oj c . However, we can still apply the Karush-Kuhn-Tucker (KKT) conditions to this 
problem to obtain the necessary conditions for optimality [20]. The Lagrangian of this problem is given 

'Different types of ideal filters, i.e., high-pass, band-pass, band-stop filters, yield essentially the same result as the ideal 
low-pass filters when the bandwidth of passband is the same. 
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by 

C = ^f> g2 (l+(fe^^^ (43) 

+ X (^J S > '< J " :>ll ~ (44) 

+ j " S 2 {a 3 S^") + v 2 )du - P r ) , (45) 

where A and v are non-negative dual variables. Due to the complementary slackness, either A = or 
the source uses full power, i.e., ^ S Xs {e^ w )duj = P s ; and either v = or the relay uses full power. 
Suppose that the source does not use full power. Then, the source can increase the PSD over the frequency 
band [oj c , tt] without changing the PSD over [0, u c ). Then, the power constraint (41) at the relay is not 
affected and the rate in (39) increases. Thus, the source should use full power to yield maximum rate, 
and we have A > due to the complementary slackness. However, the relay may or may not use full 
power depending on the channel condition. The source PSD solution to the KKT conditions is given by a 
modified water-filling method: 

( 1 (b 2 S 2 + l)a 2 \ + | | 

S x (e> u ) = { V(21n2)(A+i^) J ' M < Wc ' (46) 

a 2 ) , uj r < \w\ < 7T, 



(2 In 2) A 

where 7+ := max{0, 7}. The optimal source PSD is given by the difference between the 'water level' 
and the overall effective 'noise level'. The water level is a function of two dual variables A and v, and 
we may have two water levels if v > 0, i.e., the relay also uses full power. The effective noise level 77 is 
defined as the inverse of the overall CNR density, and it is given by r] pass = /fq^gp cr 2 and rj st0 p = cr 2 
for [0, u c ) and [lo c , tt], respectively. From here on, we will consider only the case r) paS s < cr 2 . (Otherwise, 
it is better not to use the relay (5 = 0) since a 2 over [0, tt] is the noise level without the relay.) To obtain 
optimal 5, we need to consider both (39) and (41). By differentiating rj pass (5) = w.r.t. 5 and 



setting the derivative to zero, we obtain S = % for the minimum noise level tj* 3 = ( 1 + a 2 ) a 2 / ( 1 + a 2 ) 2 . 



(1+abS)' 2 

Ipass 



However, we have the relay power constraint (41), yielding 6 < y P T j Jq C a 2 S Xg (&> u )cku + ^cr 2 ) = 
\J P r / {a 2 Ppass + (w c /7r)cj 2 ), where P pass is the portion of the source power allocated to the relay's 
passband [0, 07) and P s t op '■= P s — Ppass- Thus, the optimal 5* is given by 




Pr 



(L 2 Ppass w CJ 2 



(47) 



since r) pa ss{d) is monotonically increasing as 5 decreases from a/b to zero. With the optimal 5*, the 
multiple-tap relay filter generates a well in the noise level, as shown in Fig. 4. The effective noise level 
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Vpass ($) of this well is equal to or lower than that of the AF scheme because the effective noise level for 
the AF scheme is i]af{cL) = ^j^^p for — tt < ui < tt (see (24)) and because for the same source 
power P s and relay power P r the upperbound for 5 in (47) is larger than that for d in (25). (Note that 
Ppass < Ps and io c < tt.) Thus, when P r is small or b << o such that (s.t.) d* < a/b (consequently, 
Vpass(8) < f]AF{d)), and P s is also small enough to be confined within the passband well, the ideal low- 
pass filtering outperforms the AF scheme. (The amount of gain will be evaluated numerically shortly.) 

The structure of solution S Xs {e^) to the KKT conditions can be classified into four types depending 
on whether the source power is allocated to the relay's stopband [lo c , tt] or not, and whether the relay uses 
full power or not (equivalently, 5* / a/b or not). Fig. 4 shows the solution types. A solution of Type 1 



max J£ log 2 (l + T £+^ . £ Pj)ass ) + *=e* i og2 (! + £•- 



2— P. 



stop 



s.t. P n 



YES 



Type 1-1 




(A, !/) = (+,<>) 

{Ppass •> Pstop) — (~f~} 0) 



(A, i/) = (+,+) 

{Ppass i Pstop*) — (H~) 0) 



(A,i/) = (+,0) (A, !/) = (+,+) 

{Ppass i Pstop) — (H - ) ~t~) {Ppass i Pstop) — ("h? H~) 



— ?]pas 



J ijpass 



/2 



J Tjpass 



r 2 



Vpa 
OJc 



h 

2 



Fig. 4 



Different types of the solution to the KKT conditions (li = (21n2)(A 1 +iya2t52) and l 2 = (21n2)A , 



J. ) 

In 91 X / 



occurs when all source power is allocated to the passband of the relay filter. We can further distinguish 
Type 1 depending on the power use of the relay. When the relay uses full power, both Lagrange dual 
variables A and v are positive, and there exist two water levels, although no water or power is allocated 
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in the stop band, i.e., I2 < a 2 (Type 1-2). Otherwise, we have Type 1-1 in which only A is non-zero. 
A solution of Type 2 occurs when the relay does not use full power and the source power is allocated 
over the entire band; there is one water level common to both the passband and stopband. A Type 3 
solution occurs when the relay uses full power and the source power is distributed over the entire band; 
both Lagrange dual variables are positive and the water levels at the passband and stopband are different. 
Different types of solutions occur for different combinations of parameters, P s , P r ,a 2 ,a, b, and ui c . 




Transmit power (P =P ) Transmit power (P =P ) Transmit power (P =P ) 

(a) (b) (c) 

Fig. 5 

Rate performance: Ideal low-pass relay filtering versus AF (P s = P r , a 2 = 1, a = 1 and b = 2): 
(A) FIXED UJ c (at low P s ), (b) FIXED U c , and (C) OPTIMIZED UJ c 




To evaluate the performance of the ideal low-pass filtering at the relay, the optimization problem (39 - 
42) was solved numerically using a commercial tool with S and P paS s as variables for given P s , P r , a 2 , 
a, b and uj c . First, we fixed a = 1, b = 2 and a 2 = 1 (the same as for Fig. 3 (a) and (b)), and swept 
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P s = P r for each ui c G {0.37T, 0.67r}, and compared the performance of the ideal low-pass filtering with 
that of the AF scheme. The performance is shown in Fig. 5 (a) and (b). (The same rate curves are plotted 
with two different x-axis ranges.) It is seen in Fig. 5 (a) that indeed there is a gain over the AF scheme at 
the very low source power values. This is because the low-pass relay filter intentionally generates a well 
in the noise level in the passband and because all source power is allocated into this well at the very low 
source power values, as explained already. In this way, the signal is transmitted through a narrowband 
channel that has a lower noise level than that of the AF scheme. As P s increases, however, the source 
power spills over the relay's stopband. At high SNR, the effect of the intentional noise level shaping is 
negligible (the effect of noise itself becomes negligible), and the uniform source power distribution over 
the entire frequency band is optimal at high SNR. Consequently, filtering out the signal existing at the 
stopband [u c , it] at the relay is detrimental to the performance at high SNR, as shown in Fig. 5 (a) and (b). 
So, we optimized the cut-off frequency uj c together with 5 and P paS s, and the result is shown in Fig. 5 (c). 
It is seen that the optimized solution shows the performance gain over the AF scheme at low SNR and 
eventually converges to the AF scheme as P s increases. In this case of a = 1 and 6 = 2, the performance 
gain of the LPF relay over the AF scheme is not significant since the AF scheme performs well already, 
as shown in Fig. 3 (a) and (b) [6, Proposition 9]. Thus, we tried another case of a = b = 2 and a 2 = 1 
in which the AF scheme has a noticeable loss from the cut-set bound, as shown in Fig. 3 (c) and (d). 
Fig. 6 shows the achievable rate in this case. The figure shows similar behavior to the case of a = 1 and 
6 = 2. One important fact here is that the performance gain of the ideal low-pass (multiple-tap) relay 
filtering (even breaking the causality constraint) over the AF scheme is not so significant in both of the 
cases. Thus, it seems that, practically, there is no need for complicated linear (time-invariant) filtering 
at the relay beyond the AF scheme in flat-fading channels. However, the insignificant gain is only for 
flat-fading channels, and this is not the case in realistic ISI channels. In the next section, we will tackle 
the LTI relay problem (17, 19, 20, 21) in ISI channels, propose a practical method to solve this problem, 
and show that the LTI relay filtering yields a noticeable gain over the AF scheme in ISI channels. 

V. Joint Source and Relay Filter Design in ISI Channels 

The LTI relay problem (17, 19, 20, 21) in ISI channels is basically a constrained optimization problem. 
Our approach to this problem is based on a powerful tool of the projected subgradient method developed 
by Polyak [14] and Yamada et al. [15, 16]. In the next subsection, we will briefly introduce their results 
that are relevant to our problem, and then move on to the original LTI relay problem. 
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A. Background: The Adaptive Projected Subgradient Method (APSM) [ 15, 16] 

The projected subgradient method was initially proposed by Polyak [14] to solve cost minimization 
problems over convex constraint sets. Recently, Yamada et al. fully generalized the method even to 
adaptive situations with time-varying cost functions, and proved strong convergence of the method [15, 
16]. 

Definition 2 (Subdifferential) Let (p ■ H — > R be a continuous convex function from a Hilbert space % 
to the set R of real numbers. Then, the subdifferential of <f> at u is defined as the set of all the subgradients 
of 4> at u: 

d<j>(u) := {g G H\ < v- u,g > + 0(u) < 0(v), Vv G H}, (48) 
where < •, • > is the inner product between two vectors. 

Here, < v — u, g > + 4>(u) = r is a hyperplane with coordinates (v, r) passing (u, 4>(u)) in the 
space % x R. Thus, (48) implies that the hyperplane < v — u, g > + 0(u) = r is below the cost 
surface 0(v) = r. It is known that d<p(u) is nonempty, and by definition we have G d(j)(u) 4^ </>(u) = 
min ve % </>(v). In the differentiable case, the gradient is a unique subgradient. 

Definition 3 (Metric projection) The metric projection Pk(u) of a point u onto a closed convex set K 
is the closest point of u in K, i.e., ||Pr-(u) — u|| < ||v — u||, V v G K. 

Definition 4 (Subgradient projection) Suppose that a continuous convex function (f> : % — > R satisfies 
lev<o^ := {v G T-L\<p{v) < 0} ^ 0, i.e., the zero-level set is not empty. Let (j)' : % — > % be a selection 
of the subdifferential d<j>, i.e., 4>'(u) G d<p(u), Vu G T~L. Then, a mapping T sp u\ :7i—t'H defined by 

u, if 0(u) < 0, 

is called a subgradient projection relative to (p. 

Here, the normalization of the subgradient by factor pTT^jrp is crucial to determining the step size pa- 
rameter. Basically, the subgradient projection moves the current point u to its metric projection onto 
the intersection of two hyperplanes r = and < v — u, g > + <fi(u) = in the space of % x R 
with coordinates (v, r). (See [21, Fig. 16].) It is known that the subgradient projection always moves the 
original point closer to every point in the zero-level set, which is known as the Fejer monotone property 
with respect to the zero-level set [16, 22]. It is also known that the subgradient projection belongs to the 
class of firmly quasi-nonexpansive mappings [16,22]. Now, we state the adaptive projected subgradient 
theorem by Slavakis, Yamada and Ogura. 
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Theorem 5 (APSM, [15, 16]) Let <p n : % — > [0,oo) (Vn G N) be a sequence of continuous convex 
functions and K C Ti a nonempty closed convex set. For an arbitrarily given uo G K, the sequence 
{u n } ng N C K generated by the adaptive projected subgradient method: 

#.(«») 1 1 

u n , otherwise, 



u n+l — N PUJ 



where </>^(u n ) G <9(/) n (u n ) and < /i n < 2, satisfies the following. 
(Monotome approximation) Suppose 

u„ i n n := {u G K\<f> n (u) = <f>* n := inf (f> n (v)} + 0. (51) 

VGA' 

Then, by using V/i n G (0,2(1 - ^ n) )),we have ||u re+ i - u^|| < ||u„ - u^|| for all u (n) G U n . 
(Asymptotic optimality) Suppose 

3N G N s.t. <p* n = Vn > N and := n n >jv ^n ¥= 0- ( 52 ) 

Then, {u n } n6 N is bounded. Moreover, if we specially use V/x n G [ei,2 — 62] C (0,2), we have 
lim n ^oo (fi n (u n ) = provided that {^(u n )} neN is bounded. 

(Strong convergence) Assume (52) and SI has some relative interior w.r.t. a hyperplane II (c %); i.e., 
there exist u 6 Ilnfl and 3 e > satisfying {u G II s.t. ||u — u|| < e} C Q. Then, by using 
V?? n G [ei, 2 — €2] C (0, 2), {u n } ne pj converges strongly to some u G K, i.e., lim^oo ||u n — u|| = 0. 
Moreover, lim^oo </> n (u) = if {^> / „(u n )} n6 pj is bounded and if there exists bounded {^'„(u)} n6 N 
where <f>' n (u) G d(p n (u), Vn G N. 

Note that the update rule in (50) is a composite projection composed of the subgradient-based projec- 
tion relative to (j) n and the metric projection onto K. In fact, the first projection is not the subgradient 
projection exactly since the step size parameter fi n does not need to be one. This projection can be 
rewritten as 

^n( u n) : = u n - fa lllf""! ^^'') = K 1 ~ Vn)I + ^nT sp{<j>n )](u n ), (53) 
||<Pn( u nJ|| 

where / is the identity mapping. When p, n = 1, T^ n = T sp u n y When /x n < 1, the projection moves the 
original point towards the exact subgradient projection point, but not fully. When fi n > 1, on the other 
hand, the projection moves the original point beyond the exact subgradient projection point. For any fi n G 
(0, 2), T Mn is called a /i„-averaged quasi-nonexpansive mapping, and the properties of this mapping play 
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an important role in the proof of Theorem 5. The major difference of the subgradient projection method 
from the simple gradient method without normalization is that we know the exact range of the step size 
parameter fi n for convergence because of the proper normalization of the subgradient in the subgradient 
projection. (See (49) and (50).) The above result is general and can be applied to many constrained 
optimization problems. Consider the simple case with a fixed cost function (i.e., <fi n (u) = </>(u) V n) of 
which minimum is achieved in K. Let 0* = min ue if </>(u). Then, we can define <fi(u) := </>(u) — <fi*. 
Then, the condition (52) is trivially satisfied, and we have lim^^oo (p(u n ) = 0, i.e., linin^oo (p(u n ) = <p*. 
In this case, the result reduces to the Polyak's projected subgradient method for which he showed the weak 
convergence [14, Theorem 1]. (For the proof of Theorem 5, refer to [16]. For the introduction of related 
projections, see [21].) 

B. Joint Source and Relay Filter Design via the Projected Subgradient Method 

Even under the LTI formulation (17, 19, 20, 21), the search space for source and relay filters T(z) 
and H(z) has countably infinite dimensions since both T(z) and H{z) can be IIR filters. To avoid the 
difficulties in searching in an infinite dimensional space and in imposing the stability condition, we restrict 
ourselves to the case that both T(z) and H(z) have FIRs as in most practical filtering applications. Then, 
the stability and causality constraints are automatically satisfied. In this case, the source and relay filter 
responses are respectively given by 

T{z) = to + hz' 1 + --- + t Ls -iz- Ls+1 , and (54) 
H(z) = h + h lZ ~ l + --- + h Lr . lZ ~ u+1 , (55) 

where L s and L r are the orders of source and relay filters, respectively. Although the channel responses 
H sr (z), H rc [(z) and H s d(z) do not need to have finite durations, we also assume that these channel 
responses have finite duration L for simplicity. For notational convenience, we define the following: 

w m (w) := [l ) e*V 2w > ... ,e^ m -^] T , (56) 
t := [in, ti,"- ,*L s -i] T £» LsXl , (57) 

h := [ho,hi,'~ M r ~i] T eR irXl , (58) 

u := [t T ,h T ] T Gl^ 1 ^ 1 , (59) 
Kb ■= [h a b[0},hab[M, ■■■ ,h a b[L - 1]] T for (a, 6) = (s,r), (r,d) and (s,d). (60) 
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Then, we have 

T(e^)=w£>)t, tf(e^) = w£»h, H ab (e?») = wf (u;)h a6 . (61) 

Here, we assume that t and h are real vectors for simplicity, but the extension to the complex case is 
straightforward. Now from (17) we define the cost function as 

1 f !, (, , + H, r (e'")H(e'»;h)H rd (z'")\ 2 m ,„ ^ a 



« u > - -SF L 5 log2 1 1 + "^WWT! |T(e " ; 4)1 J <62) 

where the respective dependence of H(e ,u} ) and T(e JaJ ) on h and t is explicitly shown. The gradient 
<fi'(u) at u can be obtained after some computation as 

i r i i / \ 

where A(w) = CNR(e^;h) • t T w Ls w£t, = CNR(e JW ;h) • (w Ls wf s + w^w^Jt, d(u) = 

t T w Ls w£t /£>(e*";h), C 2 (w) = V h A7(e^;h), C 3 (u) = V h D(e^;h) and C a (cj) = CNR(e^;h); 
and N(e? u ;h) and L)(e- JCJ ;h) are the numerator and denominator of CNR(e : ' w ;h), respectively. The 
source power constraint (19) is given in terms of t and h by 

±- ^ \T{en\ 2 du = £ tf = ||t|| 2 < P s (64) 
n ^ 7r 1=0 

by Parseval's theorem, and the constraint set for t determined by (64) is simply a ball in t, denoted by 
B t (P s ), with no constraint on h. Next, consider the power constraint at the relay. The relay power 
constraint (20) is expressed in terms of t and h as 

^ J* \H(e? u )\ 2 (\H„(e? u )\ 2 ■ |T(e^)| 2 + a 2 )doo (65) 

= _L f h r w L >)w£»h • (|if sr (e^)| 2 t T w Ls ( W )wf s (a;)t + a 2 ) du;, 

/7T J~7T V v ' 

=:g(u;t) 

= hT (^/ ^^(wJwjHdwjh < P r . (66) 

Note that the constraint set for (t, h) determined by (66) is bi-convex, i.e., convex in each of t and h but 
not jointly convex. The constraint set for h for a given t is an ellipsoid, as shown in (66). Let us denote 
this ellipsoid by £h(t). The above inequality can also be written as 

t T |^|^ r .(e^)| 2 h T w Lr ( W )wf r (a;)h)w Ls ( W )wf s (^)da;) t 

+ ^y\ T w L >)wf r Hh^ < P r , (67) 
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and, therefore, the constraint set for t for a given h is also an ellipsoid. 
Now, the problem under the LTI FIR formulation is given by 

min0(t,h) (68) 

t,h 

such that 

(t,h) G Ci and (t,h) e C 2 , (69) 

where </>(t,h) is given by (62), and Ci and C2 are the constraint sets determined by (64) and (66), 
respectively. The situation is depicted in Fig. 7. Even though </>(t, h) is not jointly convex in t and h, 
we can still apply the projected subgradient method in Theorem 5. Suppose now that K := C\ n C2 is 
convex. Then, applying the projected subgradient method (50) with uo 6 K leads to the convergence to 
a local minimum at least. Due to the nonconvexity of C2, however, K is not a convex set. To circumvent 
this issue, let us investigate the structure of the two constraint sets further. Fig. 8 shows the two constraint 
sets, C\ and C2, and their intersection K in the case of t = (to, t\) and h = ho- C\ and C2 are the red 
cylinder and the blue shape that looks like a mountain on one side and is symmetric about the ho = 
plane in the left side of Fig. 8, respectively. In this case, the ellipsoid £h(t') is a line segment passing 
through (t',0) ending at the upper and the lower surfaces of the blue shape, as shown in Fig. 8. As 
mentioned already, for a given h the constraint set for t is an ellipsoid. (See (67).) So, we have an 
ellipsoidal cut of the blue shape by a plane perpendicular to the /io-axis. Thus, C\ n C2 is given by the 
shape in the right side of Fig. 8, which looks like a cylinder with bulges on the top and the bottom. If we 
are willing to sacrifice the bulging regions from the feasibility set, we can construct a convex set K C K 
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ho 




Fig. 8 

The structure of two constraint sets in case of t = (t , ti) and h = ho 

as the Cartesian product of B t {P s ) and Hh(-P s , P r ), i.e., 

K:=B t {P s ) xE h (P s ,P r ), (70) 

where 

s h (^, J Pr)= n eh(t)= n £h(t). 

teB t (P s ) tesurface(B t (-P s )) 
The second equality in (71) is because £h(t') C £h(7*') f° r 7 £ [0, 1]. This can easily be verified from 
(65). For -ft', \T(e ju) )\ 2 will simply scale down to -y 2 \T(e ju; )\ 2 , and \H{e ju) )\ 2 satisfying the inequality 
with \T{e^)\ 2 will satisfy the inequality with 7 2 |T(e- Ja; )| 2 . Here, r,\ l {P s ,P r ) is convex since it is the 
intersection of ellipsoids. Also, K is convex due to its Cartesian product structure. In the considered case 
of (to, t±, ho), K is exactly the red cylinder in the right side of Fig. 8. Based on the above discussion, we 
now provide our first algorithm for the joint source and relay filter design. 

Algorithm 1: • Initialization: Set uo properly, e.g., 

u = [t^,O r ] T (72) 
with t = [J^[l, 1, ■ ■ ■ , 1] T . (It is easy to verify that this u G K C K.) 
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• Update: 

_ ( p k ( Un - ^^^K)) , if / o, 

u n+l — \ yi j) 

u n , otherwise, 

where 4>(u) and 4>'(u) are given by (62) and (63), respectively, < fi n < 2 and u n = [t^, h^] T . 

• Stopping rule: Stop the update either when u n does not change further or when the number of iterations 
exceeds a preset maximum. 

Due to the Cartesian product structure of K in (70), the metric projection P% can be implemented by 
the composition of two separate projections, one in t projecting onto B t {P s ) and then the other in h 
projecting onto Hh(-P s , Pt), i.e., 

PK = Pz h (P s ,P t )°PB t (P a ), (74) 
where the projection onto a ball is simply given by 

t, if lltll 2 < P s , 

p w*H J (75) 

[ otherwise, 

and -P= h (p Si p t ) can be approximated by using the successive projection method [22]. That is, we uni- 
formly partition the surface of the L s -dimensional ball with radius \/P~ s into M subsets, and select the 
center t m of subset m. Then, Eh(P s , P t ) ~ n^ =1 £h(t m ) and 

PE h (Ps,P t )( h ) « ^ h (tM) ° • • • ° n h ( tl )( h )' ( 76 ) 

where the projection onto an ellipsoid can easily be implemented by a known method like the one in [23]. 
If 0'(u) ^ for all u G if, then Algorithm 1 will find the point that yields the maximum rate within K. 
Otherwise, Algorithm 1 finds a local optimum attracting uo, and convergence is guaranteed. However, the 
complexity of the projection P Sh r Ps p t \(h.) is prohibitive even for small values of L s . Thus, we propose 
a simplified algorithm to eliminate this difficulty in the following. 

Algorithm 2: • Initialization of uq. 

• Update: 

P ?h(P Bt( P s) (tn))( h n) ° PB t (P s ){tn) o (u n - , if 0nK) + 0, 

u„, otherwise. 



u n +i 



(77) 



Stopping rule: The same as that in Algorithm 1. 
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In the initialization step, we can consider other initial points than the example in Algorithm 1. For 
example, 

t = h/^,0,--- ,0] T and h = P fh(to) ([l,--- ,1] T ). (78) 

With such an initial point, we can start the algorithm from the full power use at the relay. In Algorithm 
2, the subgradient projection is the same as that in Algorithm 1, but the projection to K is now different. 
Here, we first project only the t n coordinates to the ball B t (P s ), and then project the h n coordinates onto 
the h-ellipsoid £h(^s t (p s )(tn)) corresponding to P Bt ^ Ps ^(t n ) given by (66); several projection examples 
in the case of L s = 2 and L r = 1 are shown in the right side of Fig. 8. In this way, the metric projection 
to K is approximated and highly simplified, and u n £ K for all n. The convergence of Algorithm 2 is not 
guaranteed theoretically, but there is no loss in the feasibility set in this case. It is observed numerically 
that Algorithm 2 is stable and works well. This is because K seems almost convex as seen in the example 
in Fig. 8, and the proposed two-step projection onto K approximates the metric projection onto K well 
due to the almost cylindrical structure of K except the top and bottom. Under our formulation, it is 
straightforward to impose the strict causality constraint on the relay filter H (z) that captures the possible 
delay at the relay for digital processing; simply remove ho in (55) and follow the same procedure as the 
causal case. 

VI. Numerical Results 

To evaluate the performance of the proposed joint filter design method presented in Section V, we 
ran extensive simulations under various channel conditions. The channel order L of all the propagation 
channels H sr {z), H rd {z) and H sd {z) was selected as L = 5, and each channel tap coefficient was 
generated independently and identically-distributedly (i.i.d.) according to a Rayleigh distribution with a 
different variance for a different channel, i.e., 

h sd [l] Li ~ N{Q,<i 2 sd ), h sr [l] "~N(0,a 2 sr ) and h rd [l] ^ Af(0, a 2 rd ) (79) 

for I = 0, 1, • • • , L — 1. Throughout the simulations, we fixed a 2 = 1 and used L s = 30 and L r = 20 
for the orders of the source and relay filters, respectively, to allow enough freedom for the two filters 
considering L = 5. We used Algorithm 2 with the initial point** (78) for the simulations, and the step size 
for the algorithm was adaptively changed according to /i„ = 1/ \fn as the number of iterations increases. 
The stopping condition for the algorithm was either that the square of the normalized difference of two 
** We observed that Algorithm 2 with the initialization (72) did not make a noticeable difference. 
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successive updates is less than 10~ 5 (i.e., ||u n+ i — u n || 2 /||u n || 2 < 10 -5 ) or that the number of iterations 

exceeds 1000. For the numerical integration to compute the quantities (62, 63, 66), we used a 512-point 

Gaussian quadrature method over [— 7r, tt]. Since it is not clear how to optimally apply the known coding 

strategies such as DF and CF to the case of ISI channels, we do not consider these schemes in this section, 

and thus we use the AF scheme as the performance reference. 

First, Fig. 9 shows the frequency responses related to the proposed method. The channel coefficients 

for this figure were generated randomly according to a 2 sr = a 2 rd = a 2 sd = 1 and given by H sr = 1.8833 + 

0.32542- 1 -0.0952z- 2 +0.0312z 3 -0.6138z" 4 , H rd = -0. 0728+1. 31482" 1 +0.9783z~ 2 + l. 7221z 3 - 

0.4123z~ 4 and H sd = -0.8864 - 1.8402Z" 1 - 1.6282z" 2 - 1.1738z 3 - 0.4154z" 4 . P s = P r = 1. The 

flattt dashed line represents the input spectrum of the conventional AF scheme, and the dotted line is the 

inverse of the CNR density, i.e., the effect noise level, in the case that the AF scheme is used at the relay. 

It is seen that there are two peaks in the noise level for the AF scheme around the normalized frequency 

values of 0.6 and 1. Note that the peak around 0.6 is very high and its width is also large, whereas the 

peak around 1 is mild. Now the red solid line marked with + denotes the effective noise level generated 

by the relay filter designed by the proposed method. We recognize that the designed relay filter is smart. 

One could imagine that a reasonable relay filter would suppress the frequency band having a bad overall 

response with the AF scheme and enhance the frequency band having a good overall response with the AF 

scheme so that the water-filling by the source filter might be enhanced. This is exactly the case with the 

mild noise peak around 1 ; the designed relay filter reinforces the mild noise peak. However, the designed 

relay filter suppresses the severe noise peak with a large width around 0.6 instead of reinforcing it. This 

is because the width of the severe peak is large and, thus, the loss of this frequency band will reduce the 

transmission rate. To this optimally shaped noise level, the water-filling-type source power allocation is 

performed by the designed source filter satisfying both the source and relay power constraints jointly with 

the designed relay filter; for the frequency bands around the two noise peaks the power is not allocated, 

as seen in the figure. (See the green solid line without any markers.) Similar behaviors are observed in 

other settings even though they are not shown in this paper. 

^The reason why the flat input spectrum is used for the conventional AF scheme is that even for a given relay filter, e.g., the 
AF filter, the optimal power allocation to maximize the transmission rate (17) not only with the source power constraint (19) 
but also with the relay power constraint (20) in ISI channels is not a trivial problem; this is not a simple and known water-filling 
problem with a single total power constraint because of the term \H sr (e^ u ) 2 in (20) unless H sr {e' UJ ) is constant over [— 7r, tt]. 
As far as we know, this problem was not handled before, and the two-step projection method in this paper provides one way to 
accomplish this in general ISI channels. 
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Fig. 9 

The frequency responses related to the proposed method 



Next, let us see the performance of the proposed method in various channel settings. Fig. 10 shows the 
transmission rate achieved by the proposed method averaged over 100 independent channel realizations 
for each of the different channel gain settings when P s = P r . We considered both causal and strictly 
causal filtering for the proposed method. For the strictly causal filtering, we eliminated the to and ho 
terms in the formulation and kept the same filter orders L s and L r as those in the causal filtering case. 
We considered two sets of channel conditions; in one set the direct S-D link has a reasonable strength 
compared with the S-R and R-D links (Fig. 10 (a), (b), (c), (d)), and in the other set the direct link 
is weak compared with the S-R and R-D links (Fig. 10 (e) and (f)). (We did not consider the case 
that the direct link is stronger than the S-R and R-D links since it is unnecessary to use the relay in 
this case.) It is seen in Fig. 10 that the gain obtained by the proposed joint source and relay filtering 
over the AF scheme is noticeable. In particular, Fig. 10 (c) shows the performance when P s = P r , 
a 2 = 1, a 2 d = 1, a 2 r = 1 and a 2 d = 4, which is equivalent to the flat-fading case of a = 1 and 
b = 2 in Fig. 3 (a) and (b) and Fig. 5. Compared with the flat-fading case shown in Fig. 5 (c), the 
gain by the proposed joint design over the AF scheme in ISI channels is significant. It is also seen 
that the loss caused by the strict causality of filtering is not significant in the proposed scheme. This 
is because removing one tap out of 30 or 20 taps does not reduce the degree of freedom of the filters 
much. Fig. 10 (e) and (f) show the case that the direct link is weak compared with the S-R and R-D 
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- Proposed - causal (30,20) 

- Proposed - strictly causal (30,20) 
-AF(1,1) 



10"' 10" 
Transmit power (P =P ) 




10 10" 
Transmit power (P =P ) 
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(b) 
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Fig. 10 

Transmission rate averaged over channel realizations (P s 
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1, a 2 sr = 1, a 2 rd = 10, (E) a 2 sd = 1/4, 
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links. Similar results are seen in the figure. Fig. 10 (f) shows the smallest gain among the six cases. 
This is explained as follows. As seen in (17, 18), the gain of the joint filtering depends on various 
factors. The worst situation for the performance of the joint filtering is the case that H sd {e^) « and 
Hrd(e jLU ) » such that \H rd (e ju )H(e juJ )\ » 1. In this case, the CNR density in (18) is approximated 
by CNR(eP u ) « \H sr (e^)H(e^)H rd (e^)\ 2 /(a 2 \H rd (e^)H(e^)\ 2 ) = \H sr (e^)\ 2 /a 2 ; and thus the 
impact of the relay filter disappears and only the optimal power allocation by the source filter under the 
two power constraints is effective. Even in this case, a non-negligible gain is observed. Although the 
results in the cases of P s = 2P r and 2P S = P r are not shown in this paper, we observed similar results to 
the case of P s = P r . 

VII. Conclusions 

We have considered the linear Gaussian relay problem. By adopting the LTI relay filtering and realiz- 
able input spectra, we have converted the problem to a joint design problem of source and relay filters. 
We have investigated the performance of this joint LTI filtering in flat-fading relay channels, and have 
shown the optimality of the AF scheme within the class of one-tap filters. In general ISI relay channels, 
we have developed a practical method for the joint filter design to maximize the transmission rate based 
on the projected subgradient method under the LTI FIR filtering framework, and have shown numerically 
that the gain of the proposed design is noticeable compared with the AF scheme in ISI relay channels. 
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